Exploration spatio-temporelle d’objets géographiques ponctuels
L’exemple du patrimoine architectural toulousain
Exploration spatio-temporelle d’objets géographiques ponctuels
Introduction
Cette fiche rend compte d’une série de traitements permettant d’explorer la base de données MERIMEE qui consigne l’ensemble des bâtiments classés aux monuments historiques. Nous prendrons l’exemple des fiches “Mérimée” du patrimoine architectural toulousain disponibles sur le site data.gouv.fr.
1 Préparation de l’analyse
1.1 Librairies utilisées
tidyverseregroupe plusieurs librairies pour la manipulation de données, dontdplyrsfpermet de réaliser la manipulation et le traitement de données géographiquesrasterpermet la manipulation d’information stockée en mode raster (matrice de pixel)pracmapermet de réaliser des lissages par moyennes mobiles de différents typesforecastoffre un ensemble de fonctions d’ajustement de séries temporellesggplot2pour la représentation graphique. C’est une librairie de référence avec sa popre syntaxeggspatialpour représenter des données spatiales. Syntaxeggplot2ggmappour faire des cartes. Syntaxeggplot2gganimate(et sa dépendancegifski) permet de réaliser des cartes animéescowplotoffre des fonctionnalités pour combiner plusieurs graphiquesrasterVispermet la visualisation de rasterleafletest une adaptation de la librairie javascript du même nom, permet de réaliser des cartes interactivesosmdatapermet de télécharger et d’utiliser des données OpenStreetMapRColorBrewerpropose l’utilisation et la création de palette de couleursplotlypermet la réalisation de représentations graphiques interactives
Création d’un vecteur comportant tous les noms des packages nécessaires :
my_packages <- c("sf",
"dplyr",
"tidyverse",
"pracma",
"forecast",
"ggplot2",
"cowplot",
"leaflet",
"raster",
"rasterVis",
"ggspatial",
"ggmap",
"osmdata",
"gganimate",
"gifski",
"RColorBrewer",
"plotly")La fonction installed.packages() liste tous les packages déjà installés. Utilisons-la pour récupérer la liste des packages que vous allez devoir installer.
Utilisons install.packages() pour installer les packages depuis le CRAN (Comprehensive R Archive Network). Le test conditionnel if() vérifie qu’il vous manque bien des packages avant de lancer l’installation.
Chargement de toutes les librairies à l’aide de la fonction lapply() qui permet d’appliquer la fonction library() à mon vecteur de noms de packages :
1.2 Les données
Télécharger les donnéesDeux couches géographiques en format shape sont mises à votre disposition :
1) communes.shp : Limites des communes de la métropole toulousaine récupérées sur data.gouv.fr
2) base-merimee.shp : Extraction de la BD Mérimée mise à disposition sur data.gouv.fr. Il s’agit d’une couche géographique d’objets ponctuels. Chaque point est un monument du patrimoine architectural toulousain précisément géolocalisé dans l’espace.
Chaque monument est caractérisé par une cinquantaine de variables. Extrait :
1.2.1 Import
Les données géographiques stockées en format shapefiles s’importent très facilement avec la librairie sf.
Ces données sont alors converties en objet sf (simple feature dataframe). Il s’agit tout simplement d’un tableau (dataframe) dans lequel chaque élement/ligne est associé à une géometrie. Ici, la colonne geometry (cf. ci-dessus) liste les coordonnées x/y des monuments localisés dans l’espace par un point.
Utilisez st_read() pour importer les données.
1.2.2 Caractéristique spatiale
La fonction summary() permet d’obtenir un résumé de l’ensemble du tableau de données ou de certaines variables d’intérêt. Ici, elle permet de décrire la variable geometry qui contient la forme et la localisation de chaque élément du tableau.
POINT epsg:4326 +proj=long...
3174 0 0 1.2.3 Caractéristique temporelle
La fonction colnames() renvoie le nom de l’ensemble des colonnes du tableau de données.
[1] "chpuser" "chpdeno" "chpgenr" "chppden" "chpvoca"
[6] "chpappl" "chpactu" "chptico" "chppart" "chprefp"
[11] "chpcoll" "chpreg" "chpdpt" "chpcom" "chpinsee"
[16] "chpploc" "chpaire" "chpcant" "chplieu" "chpadrs"
[21] "chpedif" "chprefe" "chpcada" "chpzone" "chpcoor"
[26] "chpcoorm" "chpcoorlb93" "chpcoormlb9" "chpcoorwgs8" "chpcoormwgs"
[31] "chpimpl" "chphydr" "chpscle" "chpscld" "chpdate"
[36] "chpjdat" "chpautr" "chpjatt" "chppers" "chpremp"
[41] "chpdepl" "chphist" "chpprot" "chpdpro" "chpppro"
[46] "chpapro" "chpmhpp" "chpsite" "chpinte" "chprema"
[51] "chpobs" "geometry" La variable nommée chpdate contient les informations sur la date de construction des bâtiments classés aux monuments historiques.
Appliquons les fonctions class(), head() et summary() à cette variable pour en apprendre un peu plus.
[1] "character"[1] NA "1776 ; 1866" NA "1883" "1881"
[6] NA Length Class Mode
3174 character character Cette variable contient de nombreuses valeurs manquantes et plusieurs dates sont parfois renseignées.
Pour ces raisons, R a automatiquement typé la variable en factor (catégorie) et non pas en numeric. Ainsi, la fonction summary() renvoie le nombre d’individus regroupés par catégories (levels) plutôt que de fournir un résumé de ses paramètres centraux et de dispersion.
1.3 Pré-traitements
Dans le cadre de cette exploration, nous allons supprimer les monuments sans date de construction et ne garder que la première date renseignée lorsque plusieurs périodes de construction sont proposées.
La multiplication des opérations sur un même objet peut être facilitée grâce à l’utilisation de tidyverse, un ensemble d’extensions contenues dans une librairie, qui utilise une syntaxe particulière (Sucre syntaxique) qui rend le code plus agréable à écrire comme à lire.
Dans l’exemple ci-dessous, nous modifions la variable chpdate en deux temps :
Extraction des 4 premiers caractères de la variable chpdate avec la fonction
substr()Transformation de cette chaîne de 4 caractères en nombre grâce à la fonction
as.integer()
Puis, nous supprimons les lignes pour lesquelles l’année de construction est inconnue (NA), grâce aux fonctions filter() et is.na(). Notez que l’usage du ! signifie ici “est différent de”.
Min. 1st Qu. Median Mean 3rd Qu. Max.
1453 1874 1911 1897 1950 2009 Les bâtiments toulousains classés aux monuments historiques les plus anciens datent de 1453, alors que les plus récents ont été construits en 2009, soit une période de 556 ans.
50% des bâtiments classés ont été construits avant 1911, un quart avant 1874 et un quart après 1950.
2 Exploration temporelle
2.1 Nombre de bâtiments par année
Pour observer la distribution temporelle du patrimoine classé aux monuments historiques à Toulouse, commençons par comptabiliser le nombre de bâtiments par année de construction.
On regroupe les lignes par année, tout en comptabilisant le nombre de lignes regroupées.
Merimee_year <- Merimee %>%
group_by(chpdate) %>%
summarize(nb = length(chpdate))
summary(Merimee_year) chpdate nb geometry
Min. :1453 Min. : 1.000 MULTIPOINT :167
1st Qu.:1774 1st Qu.: 1.000 POINT : 95
Median :1874 Median : 3.000 epsg:4326 : 0
Mean :1844 Mean : 4.256 +proj=long...: 0
3rd Qu.:1941 3rd Qu.: 6.000
Max. :2009 Max. :23.000 Il est intéressant de noter que le type d’objet sf est conservé malgré les importantes manipulations réalisées. Les géométries ont même été prises en compte. La variable geometry stocke désormais un ensemble de points (MULTIPOINT) décrivant la localisation de tous les batîments construits une même année.
Les fonctions de la librairie ggplot2 vont nous permettre d’en faire une représentation graphique.
ggplot()permet de construire un graphique croisant année de construction et nombre de bâtiments concernés.geom_line()trace une courbe reliant l’ensemble des points du graphique.labs()permet d’ajouter titre, sous-titre, noms des axes…
ggplot(Merimee_year, aes(chpdate, nb)) + geom_line(color = "orange") +
labs(x = "", y = "Nombre de bâtiments",
title = "Répartition des bâtiments classés aux monuments historiques par année de construction") 
2.2 Moyennes mobiles symétriques
Afin de lisser la courbe et de mieux visualiser les tendances, nous allons calculer des moyennes mobiles d’ordre k.
L’estimation de la moyenne en un point par la moyenne des valeurs qui l’entourent nécessite des dates équidistantes. Nous allons donc ajouter les années manquantes au tableau et y associer la valeur 0 (absence de bâtiments classés construits à cette date).
Pour cela, nous créons un tableau all_years contenant la colonne chp_date et recensant l’ensemble des entiers (années) compris entre la valeur minimale (année la plus ancienne) et la valeur maximale (année la plus récente) du tableau Merimee_year.
Puis, nous réalisons une Jointure entre les tableaux all_years et Merimee_year en utilisant leur variable identique chp_date et la fonction merge().
Les années manquantes ont été ajoutées, mais elles affichent la valeur NA (not available) dans la variable nb. Remplaçons ces valeurs NA par 0, puisqu’aucun batîment construit ces années-là n’a bénéficié de la certification.
Notons que ces nouvelles lignes ne comportent pas non plus de géométrie :
La librairie pracma et sa fonction movavg() permettent de calculer des moyennes mobiles à partir d’une série de valeurs. Les arguments de la fonction sont :
- x = la série statistique
- n = amplitude de la fenêtre temporelle (nombre de voisins)
- type = moyenne mobile à calculer. La valeur “s” correspond à la forme la plus simple. Elle calcule la valeur moyenne des n valeurs précédant la valeur courante en tenant compte de cette dernière.
Il s’agit donc d’un lissage linéaire puisqu’un poids identique est donné à chacune des valeurs considérées.
# Le résultat est stocké dans la variable ave_nb
Merimee_year$ave_nb <- movavg(Merimee_year$nb,
n= 5,
type="s")Dans un objectif comparatif, calculons des moyennes mobiles selon plusieurs valeurs de voisinage k (5, 10, 20 et 40) et représentons les différentes courbes obtenues côte à côte.
Pour cela nous utiliserons la boucle for() qui permet de reproduire, en boucle, une serie de traitements. A chaque nouvelle boucle, nous modifions la valeur de k et réalisons une représentation graphique qui est enregistrée dans un objet. La fonction assign() personnalise le nom de cet objet ggplot à l’aide de la valeur de k.
for (k in c(5,10,20,40)) {
# Calcul des moyennes mobiles
Merimee_year$ave_nb <- movavg(Merimee_year$nb, n= k, type=c("s"))
# Construction du graphique
Wind_plot <- ggplot(Merimee_year, aes(chpdate,ave_nb)) +
geom_line(color="orange") +
labs(x = "", y="Nombre de bâtiments", subtitle = paste0("Fenêtre temporelle = ", k, " ans"))
# Spécification du nom avec la valeur de k
assign(paste0("window_", k), Wind_plot)
}Une fois que les quatre graphiques ont été construits et enregistrés dans des objets différents, il est aisé de les afficher tous ensemble dans la fenêtre graphique grâce aux fonctions proposées par la librairie cowplot.
# Association des graphiques dans une grille
multi_plots <- plot_grid(window_5,window_10,window_20,window_40)
# Paramétrages titre/mise en page
title_plot <- ggdraw() +
draw_label("Répartition du nombre de bâtiments historiques par année de construction",
fontface = 'bold', x = 0,
hjust = 0) +
theme(plot.margin = margin(0, 0, 0, 7))
# Affichage de l'ensemble
plot_grid(title_plot, multi_plots,
ncol = 1,
rel_heights = c(0.1, 1))
2.3 Lissage exponentiel simple
Le lissage symétrique par moyennes mobiles donne un poids équivalent à l’ensemble des années contenues dans la fenêtre temporelle de dimension k. Une alternative consiste à faire décroître le poids de ces dernières en fonction de leur distance temporelle à la date d’intérêt. Il s’agit d’un lissage exponentiel.
Différents modèles de lissage exponentiel existent. La version la plus simple consiste à estimer une valeur à un temps t de la façon suivante :
\(\hat{y_t} = \alpha y_t + (1-\alpha)y_{t-1} + (1-\alpha)^2 y_{t-2} + (1-\alpha)^3 y_{t-3} + ... etc.\)
La décroissance des pondérations en remontant dans le temps est définie de façon exponentielle. Cette méthode est plus réactive aux courtes variations de tendance.
L’utilisation d’un modèle exponentiel simple semble adapté à la distribution temporelle des années de construction des bâtiments toulousains classés aux monuments historiques puisque celle-ci ne présente pas de saisonnalité ou une tendance unique.
La fonction ses() (simple exponential smoothing) de la librairie forecast permet d’ajuster un modèle exponentiel simple à une serie statistique. L’ajustement consiste à déterminer la valeur du coefficient de lissage \(\alpha\).
# Création d'un objet time series
Merimee_ts <- ts(Merimee_year$nb, start = min(Merimee_year$chpdate), end = max(Merimee_year$chpdate))
Merimee_ts_fit <- ses(Merimee_ts)
Merimee_ts_fit$modelSimple exponential smoothing
Call:
ses(y = Merimee_ts)
Smoothing parameters:
alpha = 0.292
Initial states:
l = 0.2912
sigma: 2.1595
AIC AICc BIC
4383.285 4383.329 4396.253 La valeur du paramètre \(\alpha\) est de 0.29, ce qui signifie que 29% des prédictions sont basées sur l’observation la plus récente. Les données anciennes ont un poids relativement important sur les valeurs prédites à un temps t (\(\alpha\) est plus proche de 0 que de 1), alors que les changements récents ont un impact moindre sur les valeurs prédites.
En appliquant la fonction summary() sur les résultats de la fonction ses(), on peut afficher les valeurs prédites (forecats) pour les 10 intervalles de temps (ici année) suivant la dernière date de sa série statistique (soit de 2010 à 2019).
Forecast method: Simple exponential smoothing
Model Information:
Simple exponential smoothing
Call:
ses(y = Merimee_ts)
Smoothing parameters:
alpha = 0.292
Initial states:
l = 0.2912
sigma: 2.1595
AIC AICc BIC
4383.285 4383.329 4396.253
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.007292547 2.155576 1.117794 -Inf Inf 0.8853182 -0.006887626
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
2010 1.477131 -1.290324 4.244586 -2.755326 5.709588
2011 1.477131 -1.405872 4.360134 -2.932041 5.886303
2012 1.477131 -1.516963 4.471226 -3.101941 6.056204
2013 1.477131 -1.624078 4.578341 -3.265759 6.220022
2014 1.477131 -1.727615 4.681877 -3.424105 6.378367
2015 1.477131 -1.827910 4.782172 -3.577493 6.531755
2016 1.477131 -1.925249 4.879512 -3.726361 6.680623
2017 1.477131 -2.019881 4.974143 -3.871087 6.825349
2018 1.477131 -2.112018 5.066280 -4.011998 6.966261
2019 1.477131 -2.201848 5.156110 -4.149382 7.103644autoplot() et autolayer() permettent d’afficher graphiquement la qualité de l’ajustement du modèle aux données observées.
La représentation graphique d’un modèle de lissage exponentiel présente par défaut les valeurs prédites par ce dernier. En mauve clair figure ainsi l’intervalle de confiance à 95% et en mauve foncé l’intervalle de confiance à 80 %.
Les prédictions par lissage exponentiel simple n’opèrent que pour un futur proche. Le graphique montre par ailleurs que ces valeurs sont estimées à un certain niveau (95% de chance qu’il n’y ait pas plus de 7 bâtiments construits après 2009 classés aux monuments historiques) mais de façon plate, puisque ce type de lissage ne tient compte d’aucune tendance ni saisonnalité.
3 Exploration spatiale
3.1 Carte d’inventaire interactive
La réalisation d’une carte d’inventaire interactive avec la librairie Leaflet permet de positionner chaque bâtiment dans l’espace toulousain et d’associer à chacun des points une fenêtre interactive.
Leaflet propose d’afficher le fond de carte principal d’OpenStreetMap et sa projection par défaut (web Mercator, code epsg:4326). Il est donc parfois nécessaire de transformer le système de coordonnées de référence des données géographiques.
Pour l’ajout de labels personnalisés, nous utiliserons l’identifiant du bâtiment classé (chpuser), le type de bâtiment (chpdeno), le siècle de construction (chpscle), le numéro et nom de rue (chpadrs) et le nom de commune (chpcom).
Commençons par créer la variable qui contient le texte à afficher dans la fenêtre interactive.
Merimee$label <- paste0("<b> Identifiant du bâtiment :</b> ", Merimee$chpuser, "<br>",
"<b> Type de bâtiment :</b> ", Merimee$chpdeno, "<br>",
"<b> Période de construction :</b> ", Merimee$chpscle, "<br>",
"<b> Adresse :</b> ", Merimee$chpadrs, ", ", Merimee$chpcom, "<br>")La fonction leaflet(), associée à un ensemble de fonctions complémentaires permet de générer la carte. Quelques exemples :
addTiles()pour ajouter un fond de carte
addMarkers()pour ajouter des marqueurs ponctuels de localisation
addScaleBar()pour ajouter une échelle
addMiniMap()pour ajouter une surface de localisation
Ces fonctions présentent de nombreux arguments paramétrables. On notera que l’utilisation de pipe (%>%), propre à la syntaxe tidyverse, est particulièrement utile ici pour clarifier la chaîne de traitement.
Merimee_map <- leaflet() %>%
addTiles() %>%
addMarkers(
lng = unlist(map(Merimee$geometry,1)),
lat = unlist(map(Merimee$geometry,2)),
data = Merimee,
popup = Merimee$label,
icon = list(iconUrl="https://upload.wikimedia.org/wikipedia/commons/1/10/Logo-monument-historique.png",
iconSize = 25),
clusterOptions = markerClusterOptions()) %>%
setView(lng = 1.433333, lat = 43.6, zoom = 13) %>%
addScaleBar(position = "bottomright", options = scaleBarOptions(metric = T, imperial = F)) %>%
addMiniMap()
Merimee_map3.2 Carte de densité carroyée
Il peut être intéressant de visualiser l’intensité des bâtiments classés aux monuments historiques dans la métropole toulousaine sous la forme d’une cartographie thématique. Plusieurs solutions existent : cartes de densité de points, interpolation spatiale, modèles gravitaires de position, etc.
Si l’on ne souhaite pas estimer des valeurs d’intensité de présence (ce qui nécessite de définir quels sont le voisinage, le type de distance, et la fonction d’impédence de cette dernière les plus opportuns à considérer dans le cas étudié), une solution simple consiste à subdiviser l’espace d’étude par une grille régulière (carroyage) et de comptabliser le nombre de monuments localisés dans chaque carreau.
La librairie raster nous permet de construire la grille, puis de comptabiliser les monuments.
# Construction de la grille régulière
# à partir de l'étendue géographique couverte par la localisation des bâtiments
Merimee_grid <- raster(extent(Merimee),
crs = st_crs(Merimee)$proj4string, nrow=50, ncol=50)
# Calcul du nombre de points par carreau
Merimee_count <- rasterize(Merimee, Merimee_grid, field = 1, fun = "sum")Evidemment la taille de la grille et son positionnement sur l’espace altèrent le décompte. Comme pour toute analyse exploratoire de données spatiales (ESDA en anglais), il est alors recommandé de multiplier les points de vue (c’est-à-dire d’observer les changements induits par une modification de la grille).
Le résultat est un raster (matrice de pixel).
[1] "RasterLayer"
attr(,"package")
[1] "raster"Chaque pixel (carreau) est qualifié par le nombre de bâtiments classés qu’il contient. Il s’agit d’une variable quantitative de stock, mais nous pouvons interpréter cette variable comme une densité (variable quantitative relative) car les valeurs sont calculées sur des surfaces identiques.
Réalisons une carte choroplèthe pour représenter la densité (carroyée) des monuments historiques à Toulouse.
- Import de la couche géographique des communes de la métropole toulousaine
- Sélection des limites communales de la ville de Toulouse
# Sélection par le code INSEE
Toulouse_boundary <- Communes_metropole[Communes_metropole$code_insee==31555,]- Utilisation de fonctions des librairies
rasterVisetggplot2pour la cartographie
gplot(Merimee_count) +
geom_tile(aes(fill = value)) +
scale_fill_gradient(low = "light salmon", high = "dark red", na.value="white") +
coord_equal() +
labs(title = "Densité de bâtiments classés à Toulouse",
subtitle = "Source : Fiches Mérimée du patrimoine architectural, 2019",
fill = "Nombre") +
theme_void() +
layer_spatial(Toulouse_boundary, fill=NA) +
annotation_scale(location = "bl") +
annotation_north_arrow(location = "tl")
3.3 Carte animée
La représentation cartographique de données spatio-temporelles pose de nombreux challenges, comme en témoigne la large bibliographie portant sur le sujet. L’un des outils offert au chercheur est l’animation cartographique, particulièrement populaire sur le web. Il s’agit de réaliser une carte par pas de temps, puis de les faire se succéder à l’écran chronologiquement.
Nous allons réaliser une animation par changement de localisation. Nous garderons la même symbologie tout au long de l’animation et regroupons les entités cartographiées par année de construction.
Pour réaliser une carte animée représentant les bâtiments classés par année de construction, nous utilisons la librairie ggmap.
La fonction get_map() permet de récupérer des fonds de carte (Google Maps, OpenStreetMap, Stamen Maps…). getbb() nous permet alors de préciser l’emprise géographique (bounding box) souhaitée à partir d’un nom de lieu. Nous précisons également un niveau de zoom.
# Récupération d'un fond de carte
toulouse_map <- get_map(getbb("Toulouse"), maptype = "terrain", zoom = 12)La fonction ggmap() permet de réaliser des cartographies et s’utilise avec la syntaxe ggplot2. On commence par ajouter tous les bâtiments (points).
ggmap(toulouse_map) +
geom_point(data=Merimee,
aes(x=unlist(map(Merimee$geometry,1)),
y=unlist(map(Merimee$geometry,2))),
col="tomato",size=1, shape=20)
Puis on utilise la fonction transition_time() de la libraire gganimate (syntaxe ggplot2) pour réaliser une carte animée des bâtiments par année de construction. Nous indiquons chpdate comme variable temporelle à utiliser.
La fonction shadow_mark() permet de contrôler l’affichage des entités antérieures et postérieures au pas de temps en cours et ease_aes() gère l’affichage des données durant les transitions (apparition/disparition).
ggmap(toulouse_map) +
geom_point(data=Merimee,
aes(x=unlist(map(Merimee$geometry,1)),
y=unlist(map(Merimee$geometry,2))),
col="tomato",size=1, shape=20) +
transition_time(Merimee$chpdate) +
labs(title = "Year: {frame_time}") +
shadow_mark(alpha = 0.3, size = 0.5) +
ease_aes("linear")